Spatial regression analysis for the total number of cases
load("../../data/unionedListw_d_berlinNotSubdivided.Rdata")
load("../../data/kreiseWithCovidMeldedatumWeeklyPredictors.Rdata")Explorative analysis
Scatterplot matrix
st_drop_geometry(kreiseWithCovidMeldeWeeklyCleanedPredictorsAdded) %>%
dplyr::select(incidenceTotal, Langzeitarbeitslosenquote, Hochqualifizierte, Breitbandversorgung,
Stimmenanteile.AfD, Studierende, Einpendler, Auspendler, Einpendler, Auslaenderanteil,
Laendlichkeit) %>%
ggpairs()
We see some predictors correlating with the incidence rate as well as some moderate collinearity between some of our predictors.
Maps of response and predictors
For reference I create a series of maps for the response and all predictors.
kreiseWithCovidMeldeWeeklyCleanedPredictorsAdded %>%
tm_shape() + tm_polygons(col = c("incidenceTotal"), legend.hist = TRUE, palette = "-plasma",
legend.reverse = TRUE, title = "Incidence rate for total period") + tm_layout(legend.outside = TRUE,
bg.color = "darkgrey", outer.bg.color = "lightgrey", legend.outside.size = 0.5,
attr.outside = TRUE, main.title = "Incidence rate total period") + tm_scale_bar()
mForeigner <- kreiseWithCovidMeldeWeeklyCleanedPredictorsAdded %>%
tm_shape() + tm_polygons(col = c("Auslaenderanteil"), legend.hist = TRUE, legend.reverse = TRUE,
title = "Share of foreigners", style = "pretty") + tm_layout(legend.outside = TRUE,
bg.color = "darkgrey", outer.bg.color = "lightgrey", legend.outside.size = 0.5,
attr.outside = TRUE, main.title = "Foreign population") + tm_scale_bar()
print(mForeigner)
kreiseWithCovidMeldeWeeklyCleanedPredictorsAdded %>%
tm_shape() + tm_polygons(col = c("Auspendler"), legend.hist = TRUE, legend.reverse = TRUE,
title = "Outgoing commuters as share\nof total employees", style = "kmeans") +
tm_layout(legend.outside = TRUE, bg.color = "darkgrey", outer.bg.color = "lightgrey",
legend.outside.size = 0.5, attr.outside = TRUE, main.title = "Outgoing commuters") +
tm_scale_bar()
kreiseWithCovidMeldeWeeklyCleanedPredictorsAdded %>%
tm_shape() + tm_polygons(col = c("Einpendler"), legend.hist = TRUE, legend.reverse = TRUE,
title = "Incomming commuters as share\nof total employees", style = "kmeans") +
tm_layout(legend.outside = TRUE, bg.color = "darkgrey", outer.bg.color = "lightgrey",
legend.outside.size = 0.5, attr.outside = TRUE, main.title = "Incomming commuters") +
tm_scale_bar()
mRural <- kreiseWithCovidMeldeWeeklyCleanedPredictorsAdded %>%
tm_shape() + tm_polygons(col = c("Laendlichkeit"), legend.hist = TRUE, palette = "Purples",
legend.reverse = TRUE, title = "Share of inhabitants in places\nwith less then 150 Inh/sqkm",
style = "kmeans") + tm_layout(legend.outside = TRUE, bg.color = "darkgrey", outer.bg.color = "lightgrey",
legend.outside.size = 0.5, attr.outside = TRUE, main.title = "Ruralness") + tm_scale_bar()
print(mRural)
kreiseWithCovidMeldeWeeklyCleanedPredictorsAdded %>%
tm_shape() + tm_polygons(col = c("Studierende"), legend.hist = TRUE, palette = "Purples",
legend.reverse = TRUE, title = "Students per \n1000 Inhabitants", style = "kmeans") +
tm_layout(legend.outside = TRUE, bg.color = "darkgrey", outer.bg.color = "lightgrey",
legend.outside.size = 0.5, attr.outside = TRUE, main.title = "Students") +
tm_scale_bar()
kreiseWithCovidMeldeWeeklyCleanedPredictorsAdded %>%
tm_shape() + tm_polygons(col = c("Hochqualifizierte"), legend.hist = TRUE, legend.reverse = TRUE,
title = "Highly qualified employees/ total employees", style = "kmeans") + tm_layout(legend.outside = TRUE,
bg.color = "darkgrey", outer.bg.color = "lightgrey", legend.outside.size = 0.5,
attr.outside = TRUE, main.title = "Highly qualified employees") + tm_scale_bar()
kreiseWithCovidMeldeWeeklyCleanedPredictorsAdded %>%
tm_shape() + tm_polygons(col = c("Breitbandversorgung"), legend.hist = TRUE,
palette = "Purples", legend.reverse = TRUE, title = "Share of households with internet \nconnectivity > 5 mBits/s",
style = "kmeans") + tm_layout(legend.outside = TRUE, bg.color = "darkgrey", outer.bg.color = "lightgrey",
legend.outside.size = 0.5, attr.outside = TRUE, main.title = "Highspeed internet") +
tm_scale_bar()
kreiseWithCovidMeldeWeeklyCleanedPredictorsAdded %>%
tm_shape() + tm_polygons(col = c("Langzeitarbeitslosenquote"), legend.hist = TRUE,
palette = "Oranges", legend.reverse = TRUE, title = "Unemployment rate [%]") +
tm_layout(legend.outside = TRUE, bg.color = "darkgrey", outer.bg.color = "lightgrey",
legend.outside.size = 0.5, attr.outside = TRUE, main.title = "Long term unemloyment rate") +
tm_scale_bar()
kreiseWithCovidMeldeWeeklyCleanedPredictorsAdded %>%
tm_shape() + tm_polygons(col = c("Stimmenanteile.AfD"), legend.hist = TRUE, palette = "Blues",
legend.reverse = TRUE, title = "Share of votes") + tm_layout(legend.outside = TRUE,
bg.color = "darkgrey", outer.bg.color = "lightgrey", legend.outside.size = 0.5,
attr.outside = TRUE, main.title = "Right winged party (AfD)") + tm_scale_bar()
Normal GLM
We start with a normal GLM before checking for spatial autocorrelation in the residuals. Since we have count data a Poisson GLM with an offset for the population at risk seems a natural choice.
Poisson GLM
modelGlm <- glm(sumCasesTotal ~ Stimmenanteile.AfD + Auslaenderanteil + Hochqualifizierte +
Langzeitarbeitslosenquote + Einpendler + Studierende + Laendlichkeit + offset(log(EWZ)),
data = kreiseWithCovidMeldeWeeklyCleanedPredictorsAdded, family = poisson)
summary(modelGlm)
Call:
glm(formula = sumCasesTotal ~ Stimmenanteile.AfD + Auslaenderanteil +
Hochqualifizierte + Langzeitarbeitslosenquote + Einpendler +
Studierende + Laendlichkeit + offset(log(EWZ)), family = poisson,
data = kreiseWithCovidMeldeWeeklyCleanedPredictorsAdded)
Deviance Residuals:
Min 1Q Median 3Q Max
-75.552 -13.633 -2.309 10.382 70.683
Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) -3.449e+00 4.460e-03 -773.31 <2e-16 ***
Stimmenanteile.AfD 4.122e-02 9.852e-05 418.38 <2e-16 ***
Auslaenderanteil 3.289e-02 1.226e-04 268.35 <2e-16 ***
Hochqualifizierte -1.276e-02 1.185e-04 -107.66 <2e-16 ***
Langzeitarbeitslosenquote -2.364e-02 3.704e-04 -63.81 <2e-16 ***
Einpendler -1.599e-03 4.323e-05 -36.99 <2e-16 ***
Studierende 1.509e-04 1.211e-05 12.46 <2e-16 ***
Laendlichkeit -1.327e-03 2.565e-05 -51.72 <2e-16 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
(Dispersion parameter for poisson family taken to be 1)
Null deviance: 344479 on 399 degrees of freedom
Residual deviance: 156226 on 392 degrees of freedom
AIC: 160590
Number of Fisher Scoring iterations: 4Negative binomial GLM to account for overdispersion
Since we should be suspicious with respect to overdispersion we will run a negative binomial and afterwards a quasi-poisson GLM to account for that. Since the negative binomial GLM triggers some complications when using it with the spatial eigenvector mapping we will stay with the quasi-poisson model afterwards.
modelGlmNb <- glm.nb(sumCasesTotal ~ Stimmenanteile.AfD + Auslaenderanteil + Hochqualifizierte +
Langzeitarbeitslosenquote + Einpendler + Studierende + Laendlichkeit + offset(log(EWZ)),
data = kreiseWithCovidMeldeWeeklyCleanedPredictorsAdded)
summary(modelGlmNb)
Call:
glm.nb(formula = sumCasesTotal ~ Stimmenanteile.AfD + Auslaenderanteil +
Hochqualifizierte + Langzeitarbeitslosenquote + Einpendler +
Studierende + Laendlichkeit + offset(log(EWZ)), data = kreiseWithCovidMeldeWeeklyCleanedPredictorsAdded,
init.theta = 21.58561699, link = log)
Deviance Residuals:
Min 1Q Median 3Q Max
-3.4933 -0.7228 -0.0782 0.5668 3.1286
Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) -3.5407641 0.1201151 -29.478 < 2e-16 ***
Stimmenanteile.AfD 0.0458598 0.0022967 19.968 < 2e-16 ***
Auslaenderanteil 0.0376844 0.0030124 12.510 < 2e-16 ***
Hochqualifizierte -0.0144682 0.0029854 -4.846 1.26e-06 ***
Langzeitarbeitslosenquote -0.0480101 0.0090316 -5.316 1.06e-07 ***
Einpendler -0.0012145 0.0013595 -0.893 0.372
Studierende 0.0002775 0.0002649 1.048 0.295
Laendlichkeit -0.0008254 0.0005037 -1.639 0.101
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
(Dispersion parameter for Negative Binomial(21.5856) family taken to be 1)
Null deviance: 851.01 on 399 degrees of freedom
Residual deviance: 403.29 on 392 degrees of freedom
AIC: 7156.5
Number of Fisher Scoring iterations: 1
Theta: 21.59
Std. Err.: 1.52
2 x log-likelihood: -7138.506 drop1(modelGlmNb)Single term deletions
Model:
sumCasesTotal ~ Stimmenanteile.AfD + Auslaenderanteil + Hochqualifizierte +
Langzeitarbeitslosenquote + Einpendler + Studierende + Laendlichkeit +
offset(log(EWZ))
Df Deviance AIC
<none> 403.29 7154.5
Stimmenanteile.AfD 1 812.34 7561.6
Auslaenderanteil 1 556.93 7306.1
Hochqualifizierte 1 427.05 7176.3
Langzeitarbeitslosenquote 1 431.45 7180.7
Einpendler 1 404.07 7153.3
Studierende 1 404.45 7153.7
Laendlichkeit 1 405.90 7155.1modelGlmNb <- update(modelGlmNb, ~. - Einpendler)
summary(modelGlmNb)
Call:
glm.nb(formula = sumCasesTotal ~ Stimmenanteile.AfD + Auslaenderanteil +
Hochqualifizierte + Langzeitarbeitslosenquote + Studierende +
Laendlichkeit + offset(log(EWZ)), data = kreiseWithCovidMeldeWeeklyCleanedPredictorsAdded,
init.theta = 21.54407359, link = log)
Deviance Residuals:
Min 1Q Median 3Q Max
-3.4684 -0.7223 -0.1059 0.5948 3.1251
Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) -3.6325306 0.0626247 -58.005 < 2e-16 ***
Stimmenanteile.AfD 0.0458152 0.0022989 19.929 < 2e-16 ***
Auslaenderanteil 0.0377287 0.0030153 12.513 < 2e-16 ***
Hochqualifizierte -0.0142960 0.0029816 -4.795 1.63e-06 ***
Langzeitarbeitslosenquote -0.0437307 0.0076450 -5.720 1.06e-08 ***
Studierende 0.0003129 0.0002617 1.196 0.232
Laendlichkeit -0.0008165 0.0005042 -1.619 0.105
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
(Dispersion parameter for Negative Binomial(21.5441) family taken to be 1)
Null deviance: 849.38 on 399 degrees of freedom
Residual deviance: 403.30 on 393 degrees of freedom
AIC: 7155.3
Number of Fisher Scoring iterations: 1
Theta: 21.54
Std. Err.: 1.52
2 x log-likelihood: -7139.284 modelGlmNb <- update(modelGlmNb, ~. - Studierende)
summary(modelGlmNb)
Call:
glm.nb(formula = sumCasesTotal ~ Stimmenanteile.AfD + Auslaenderanteil +
Hochqualifizierte + Langzeitarbeitslosenquote + Laendlichkeit +
offset(log(EWZ)), data = kreiseWithCovidMeldeWeeklyCleanedPredictorsAdded,
init.theta = 21.46423243, link = log)
Deviance Residuals:
Min 1Q Median 3Q Max
-3.4607 -0.7231 -0.0854 0.5758 3.1371
Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) -3.6403759 0.0624123 -58.328 < 2e-16 ***
Stimmenanteile.AfD 0.0454243 0.0022782 19.939 < 2e-16 ***
Auslaenderanteil 0.0377774 0.0030208 12.506 < 2e-16 ***
Hochqualifizierte -0.0126774 0.0026518 -4.781 1.75e-06 ***
Langzeitarbeitslosenquote -0.0424379 0.0075784 -5.600 2.15e-08 ***
Laendlichkeit -0.0008477 0.0005042 -1.681 0.0927 .
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
(Dispersion parameter for Negative Binomial(21.4642) family taken to be 1)
Null deviance: 846.24 on 399 degrees of freedom
Residual deviance: 403.31 on 394 degrees of freedom
AIC: 7154.8
Number of Fisher Scoring iterations: 1
Theta: 21.46
Std. Err.: 1.51
2 x log-likelihood: -7140.788 Quasi-Poisson GLM to account for overdispersion
modelGlmQp <- glm(sumCasesTotal ~ Stimmenanteile.AfD + Auslaenderanteil + Hochqualifizierte +
Langzeitarbeitslosenquote + Einpendler + Studierende + Laendlichkeit + offset(log(EWZ)),
data = kreiseWithCovidMeldeWeeklyCleanedPredictorsAdded, family = quasipoisson)
summary(modelGlmQp)
Call:
glm(formula = sumCasesTotal ~ Stimmenanteile.AfD + Auslaenderanteil +
Hochqualifizierte + Langzeitarbeitslosenquote + Einpendler +
Studierende + Laendlichkeit + offset(log(EWZ)), family = quasipoisson,
data = kreiseWithCovidMeldeWeeklyCleanedPredictorsAdded)
Deviance Residuals:
Min 1Q Median 3Q Max
-75.552 -13.633 -2.309 10.382 70.683
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -3.4489644 0.0884963 -38.973 < 2e-16 ***
Stimmenanteile.AfD 0.0412170 0.0019548 21.085 < 2e-16 ***
Auslaenderanteil 0.0328903 0.0024319 13.524 < 2e-16 ***
Hochqualifizierte -0.0127585 0.0023516 -5.426 1.01e-07 ***
Langzeitarbeitslosenquote -0.0236376 0.0073502 -3.216 0.00141 **
Einpendler -0.0015993 0.0008578 -1.864 0.06300 .
Studierende 0.0001509 0.0002402 0.628 0.53040
Laendlichkeit -0.0013267 0.0005090 -2.607 0.00949 **
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
(Dispersion parameter for quasipoisson family taken to be 393.7124)
Null deviance: 344479 on 399 degrees of freedom
Residual deviance: 156226 on 392 degrees of freedom
AIC: NA
Number of Fisher Scoring iterations: 4drop1(modelGlmQp, test = "F")Single term deletions
Model:
sumCasesTotal ~ Stimmenanteile.AfD + Auslaenderanteil + Hochqualifizierte +
Langzeitarbeitslosenquote + Einpendler + Studierende + Laendlichkeit +
offset(log(EWZ))
Df Deviance F value Pr(>F)
<none> 156226
Stimmenanteile.AfD 1 322415 416.9966 < 2.2e-16 ***
Auslaenderanteil 1 225877 174.7648 < 2.2e-16 ***
Hochqualifizierte 1 167833 29.1225 1.179e-07 ***
Langzeitarbeitslosenquote 1 160333 10.3051 0.001436 **
Einpendler 1 157593 3.4279 0.064854 .
Studierende 1 156381 0.3869 0.534276
Laendlichkeit 1 158920 6.7577 0.009687 **
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1modelGlmQp <- update(modelGlmQp, ~. - Studierende)
drop1(modelGlmQp, test = "F")Single term deletions
Model:
sumCasesTotal ~ Stimmenanteile.AfD + Auslaenderanteil + Hochqualifizierte +
Langzeitarbeitslosenquote + Einpendler + Laendlichkeit +
offset(log(EWZ))
Df Deviance F value Pr(>F)
<none> 156381
Stimmenanteile.AfD 1 325141 424.1105 < 2.2e-16 ***
Auslaenderanteil 1 225948 174.8296 < 2.2e-16 ***
Hochqualifizierte 1 168755 31.0982 4.57e-08 ***
Langzeitarbeitslosenquote 1 160334 9.9362 0.001745 **
Einpendler 1 157783 3.5239 0.061231 .
Laendlichkeit 1 159117 6.8761 0.009075 **
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1modelGlmQp <- update(modelGlmQp, ~. - Langzeitarbeitslosenquote)
drop1(modelGlmQp, test = "F")Single term deletions
Model:
sumCasesTotal ~ Stimmenanteile.AfD + Auslaenderanteil + Hochqualifizierte +
Einpendler + Laendlichkeit + offset(log(EWZ))
Df Deviance F value Pr(>F)
<none> 160334
Stimmenanteile.AfD 1 325149 405.0093 < 2.2e-16 ***
Auslaenderanteil 1 229927 171.0135 < 2.2e-16 ***
Hochqualifizierte 1 170097 23.9891 1.417e-06 ***
Einpendler 1 160405 0.1741 0.67669
Laendlichkeit 1 161683 3.3146 0.06943 .
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1modelGlmQp <- update(modelGlmQp, ~. - Einpendler)
drop1(modelGlmQp, test = "F")Single term deletions
Model:
sumCasesTotal ~ Stimmenanteile.AfD + Auslaenderanteil + Hochqualifizierte +
Laendlichkeit + offset(log(EWZ))
Df Deviance F value Pr(>F)
<none> 160405
Stimmenanteile.AfD 1 325303 406.0623 < 2.2e-16 ***
Auslaenderanteil 1 229951 171.2579 < 2.2e-16 ***
Hochqualifizierte 1 170969 26.0133 5.274e-07 ***
Laendlichkeit 1 161811 3.4608 0.06358 .
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Explained deviance:
1 - modelGlmNb$deviance/modelGlmNb$null.deviance[1] 0.5234146summary(modelGlmQp)
Call:
glm(formula = sumCasesTotal ~ Stimmenanteile.AfD + Auslaenderanteil +
Hochqualifizierte + Laendlichkeit + offset(log(EWZ)), family = quasipoisson,
data = kreiseWithCovidMeldeWeeklyCleanedPredictorsAdded)
Deviance Residuals:
Min 1Q Median 3Q Max
-76.612 -13.848 -2.166 11.822 72.191
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -3.6310778 0.0486591 -74.623 < 2e-16 ***
Stimmenanteile.AfD 0.0399200 0.0019120 20.879 < 2e-16 ***
Auslaenderanteil 0.0327747 0.0024445 13.408 < 2e-16 ***
Hochqualifizierte -0.0101702 0.0019856 -5.122 4.74e-07 ***
Laendlichkeit -0.0009186 0.0004915 -1.869 0.0624 .
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
(Dispersion parameter for quasipoisson family taken to be 400.4526)
Null deviance: 344479 on 399 degrees of freedom
Residual deviance: 160405 on 395 degrees of freedom
AIC: NA
Number of Fisher Scoring iterations: 4Explained deviance:
1 - modelGlmQp$deviance/modelGlmQp$null.deviance[1] 0.5343536We end up with a model of decent quality. Directions of effect seem to be aligned with expectations:
- higher share of votes for right winged party is associated with higher incidence. Presumably a proxy for the share of population opposing mask wearing and social distancing and vaccination
- higher share of foreigners is associated with higher incidence. Foreigners might not be reach by information campaigns with respect to social distancing due to language problems
- higher share of highly qualified work force is associated with lower incidence rates. For those employees it might be easier to work from home office and to avoid close contacts during work hours at office
- higher ruralness (higher share of population living in rural areas) is associated with lower incidence rates. This might be due to lower contact rates e.g. by lower share of public transport.
Checking for spatial autocorrelation in the residuals
Global Moran’s I
For regression residuals we need to use a different test as residuals are centered around zero and will sum up to zero.
lm.morantest(modelGlmNb, listw = unionedListw_d)
Global Moran I for regression residuals
data:
model: glm.nb(formula = sumCasesTotal ~ Stimmenanteile.AfD +
Auslaenderanteil + Hochqualifizierte + Langzeitarbeitslosenquote +
Laendlichkeit + offset(log(EWZ)), data =
kreiseWithCovidMeldeWeeklyCleanedPredictorsAdded, init.theta =
21.46423243, link = log)
weights: unionedListw_d
Moran I statistic standard deviate = 22.666, p-value < 2.2e-16
alternative hypothesis: greater
sample estimates:
Observed Moran I Expectation Variance
0.5078558227 -0.0003562342 0.0005027223 (moranGlmQp <- lm.morantest(modelGlmQp, listw = unionedListw_d))
Global Moran I for regression residuals
data:
model: glm(formula = sumCasesTotal ~ Stimmenanteile.AfD +
Auslaenderanteil + Hochqualifizierte + Laendlichkeit +
offset(log(EWZ)), family = quasipoisson, data =
kreiseWithCovidMeldeWeeklyCleanedPredictorsAdded)
weights: unionedListw_d
Moran I statistic standard deviate = 22.628, p-value < 2.2e-16
alternative hypothesis: greater
sample estimates:
Observed Moran I Expectation Variance
5.079319e-01 -5.494062e-07 5.038494e-04 Both model suffer from significant global spatial autocorrelation.
This implies that the usual assumption about independence of errors is violated. In turn, our standard errors might be too low, p-values too small, size (and potentially even sign) of the regression coefficients might be wrong. So we need to incorporate spatial autocorrelation in our analysis.
Plot residuals
kreiseCentroids <- st_centroid(kreiseWithCovidMeldeWeeklyCleanedPredictorsAdded,
of_largest_polygon = TRUE)Add residuals to centroids for plotting
kreiseCentroids$residGlmNb <- residuals(modelGlmNb)
kreiseCentroids$residGlmNbAbs <- abs(kreiseCentroids$residGlmNb)
kreiseCentroids$residGlmQp <- residuals(modelGlmQp)
kreiseCentroids$residGlmQpAbs <- abs(kreiseCentroids$residGlmQp)m1 <- tm_shape(kreiseWithCovidMeldeWeeklyCleanedPredictorsAdded) + tm_polygons(col = "grey") +
tm_shape(kreiseCentroids) + tm_bubbles(size = "residGlmNbAbs", col = "residGlmNb",
palette = "-RdBu", alpha = 0.9, perceptual = TRUE, scale = 0.8, border.alpha = 0.3,
title.size = "Abs residual", title.col = "Residuals", n = 3) + tm_layout(main.title = "Pearson residuals, GLM NB",
bg = "darkgrey", legend.outside = TRUE, attr.outside = TRUE) + tm_scale_bar()
m2 <- tm_shape(kreiseWithCovidMeldeWeeklyCleanedPredictorsAdded) + tm_polygons(col = "grey") +
tm_shape(kreiseCentroids) + tm_bubbles(size = "residGlmNbAbs", col = "residGlmNb",
palette = "-RdBu", alpha = 0.9, perceptual = TRUE, scale = 0.8, border.alpha = 0.3,
title.size = "Abs residual", title.col = "Residuals", n = 3) + tm_layout(main.title = "Pearson residuals, GLM QP",
bg = "darkgrey", legend.outside = TRUE, attr.outside = TRUE) + tm_scale_bar()
tmap_arrange(m1, m2)
As indicated by global Moran’s I we see that large positive and large negative residuals form some cluster. The differences between both GLM models seem neglectable.
Spatial Eigenvector Mapping
The idea behind the spatial eigenvector mapping approach is to use additional covariates that aborb the spatial autocorrelation, leading to unbiased estimators for the other predictors. The additional covariates are based on the eigenfunction decomposition of the spatial weight matrix \(W\). Eigenvectors of \(W\) represent the decompositions the spatial weight Matrix into all mutually orthogonal eigenvectors. Those with positive eigenvalues represent positive autocorrelation, whereas eigenvectors with negative eigenvalues represent negative autocorrelation. Only eigenvectors with positive eigenvalues are used for the selection.
Selection of eigenvectors
The function ME uses brute force eigenvector selection to reach a subset of such vectors to be added to the RHS of the GLM model to reduce residual autocorrelation to below the specified alpha value (defaults to 0.05). Since eigenvector selection only works on symmetric weights, the weights are made symmetric beforehand.
meQp <- spatialreg::ME(sumCasesTotal ~ Stimmenanteile.AfD + Auslaenderanteil + Hochqualifizierte +
Laendlichkeit, family = quasipoisson, data = kreiseWithCovidMeldeWeeklyCleanedPredictorsAdded,
offset = log(EWZ), listw = unionedListw_d)Refitting GLM under incorporation of the selected spatial eigenvectors:
modelSevmQp <- glm(sumCasesTotal ~ Stimmenanteile.AfD + Auslaenderanteil + Hochqualifizierte +
Laendlichkeit + fitted(meQp), family = quasipoisson, data = kreiseWithCovidMeldeWeeklyCleanedPredictorsAdded,
offset = log(EWZ))
summary(modelSevmQp)
Call:
glm(formula = sumCasesTotal ~ Stimmenanteile.AfD + Auslaenderanteil +
Hochqualifizierte + Laendlichkeit + fitted(meQp), family = quasipoisson,
data = kreiseWithCovidMeldeWeeklyCleanedPredictorsAdded,
offset = log(EWZ))
Deviance Residuals:
Min 1Q Median 3Q Max
-52.528 -10.397 -1.632 6.793 56.002
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -3.6821280 0.0447482 -82.285 < 2e-16 ***
Stimmenanteile.AfD 0.0405874 0.0018294 22.186 < 2e-16 ***
Auslaenderanteil 0.0286188 0.0022504 12.717 < 2e-16 ***
Hochqualifizierte -0.0051287 0.0018140 -2.827 0.004941 **
Laendlichkeit -0.0004370 0.0004085 -1.070 0.285383
fitted(meQp)vec2 1.2129876 0.1909764 6.352 6.03e-10 ***
fitted(meQp)vec6 -0.8132039 0.1744241 -4.662 4.32e-06 ***
fitted(meQp)vec4 -1.7032836 0.1475831 -11.541 < 2e-16 ***
fitted(meQp)vec10 0.1101595 0.1610011 0.684 0.494251
fitted(meQp)vec1 -0.7785702 0.1576694 -4.938 1.18e-06 ***
fitted(meQp)vec18 -0.7198156 0.1481683 -4.858 1.73e-06 ***
fitted(meQp)vec21 -1.4354292 0.1692912 -8.479 4.95e-16 ***
fitted(meQp)vec11 -0.3051171 0.1477679 -2.065 0.039609 *
fitted(meQp)vec19 -1.0557902 0.1680532 -6.282 9.04e-10 ***
fitted(meQp)vec5 0.5968183 0.1644262 3.630 0.000322 ***
fitted(meQp)vec8 0.5483723 0.1670307 3.283 0.001121 **
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
(Dispersion parameter for quasipoisson family taken to be 224.0348)
Null deviance: 344479 on 399 degrees of freedom
Residual deviance: 85686 on 384 degrees of freedom
AIC: NA
Number of Fisher Scoring iterations: 4The procedure added 11 spatial eigenvectors to the model. Together this leads to a more than satisfying amount of explained deviance. However, we need to keep in mind that a good share of that come from the spatial eigenvectors.
1 - modelSevmQp$deviance/modelSevmQp$null.deviance[1] 0.7512578Ruralness of the district now became insignificant so we might want to drop it from the model.
meQp <- spatialreg::ME(sumCasesTotal ~ Stimmenanteile.AfD + Auslaenderanteil + Hochqualifizierte,
family = quasipoisson, data = kreiseWithCovidMeldeWeeklyCleanedPredictorsAdded,
offset = log(EWZ), listw = unionedListw_d)Refitting GLM under incorporation of the selected spatial eigenvectors:
modelSevmQp <- glm(sumCasesTotal ~ Stimmenanteile.AfD + Auslaenderanteil + Hochqualifizierte +
fitted(meQp), family = quasipoisson, offset = log(EWZ), data = kreiseWithCovidMeldeWeeklyCleanedPredictorsAdded)
summary(modelSevmQp)
Call:
glm(formula = sumCasesTotal ~ Stimmenanteile.AfD + Auslaenderanteil +
Hochqualifizierte + fitted(meQp), family = quasipoisson,
data = kreiseWithCovidMeldeWeeklyCleanedPredictorsAdded,
offset = log(EWZ))
Deviance Residuals:
Min 1Q Median 3Q Max
-53.356 -10.321 -1.947 6.936 56.477
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -3.721873 0.036744 -101.293 < 2e-16 ***
Stimmenanteile.AfD 0.041288 0.001849 22.333 < 2e-16 ***
Auslaenderanteil 0.030151 0.002176 13.853 < 2e-16 ***
Hochqualifizierte -0.004930 0.001718 -2.870 0.00433 **
fitted(meQp)vec2 1.229067 0.193236 6.360 5.69e-10 ***
fitted(meQp)vec6 -0.791620 0.177978 -4.448 1.14e-05 ***
fitted(meQp)vec4 -1.718902 0.149249 -11.517 < 2e-16 ***
fitted(meQp)vec10 0.132188 0.163657 0.808 0.41975
fitted(meQp)vec1 -0.661805 0.146610 -4.514 8.46e-06 ***
fitted(meQp)vec18 -0.713091 0.149129 -4.782 2.48e-06 ***
fitted(meQp)vec21 -1.467720 0.171605 -8.553 2.85e-16 ***
fitted(meQp)vec11 -0.332178 0.148986 -2.230 0.02635 *
fitted(meQp)vec19 -1.001783 0.169929 -5.895 8.17e-09 ***
fitted(meQp)vec5 0.549126 0.167337 3.282 0.00113 **
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
(Dispersion parameter for quasipoisson family taken to be 229.6598)
Null deviance: 344479 on 399 degrees of freedom
Residual deviance: 88520 on 386 degrees of freedom
AIC: NA
Number of Fisher Scoring iterations: 4The eigenvectors selected did not change by excluding ruralness from the model.
Plotting selected spatial eigenvectors
Presumably we have missed a lot of other predictors that are now partially absorbed into different eigenvectors. It might be worth to plot and investigate the eigenvectors that made it into the model. Therefore, we attach the selected eigenvectors to the sf object and plot them.
summary(fitted(meQp)) vec2 vec6 vec4 vec10
Min. :-0.17702 Min. :-0.190177 Min. :-0.107585 Min. :-0.148895
1st Qu.:-0.04621 1st Qu.:-0.014842 1st Qu.:-0.030753 1st Qu.:-0.018094
Median : 0.02081 Median : 0.003256 Median :-0.008749 Median : 0.003103
Mean : 0.00000 Mean : 0.000000 Mean : 0.000000 Mean : 0.000000
3rd Qu.: 0.03945 3rd Qu.: 0.013008 3rd Qu.: 0.040463 3rd Qu.: 0.017146
Max. : 0.07715 Max. : 0.185980 Max. : 0.112574 Max. : 0.141635
vec1 vec18 vec21
Min. :-0.157693 Min. :-0.2348596 Min. :-0.1337760
1st Qu.:-0.028969 1st Qu.:-0.0206895 1st Qu.:-0.0257369
Median : 0.001129 Median : 0.0003868 Median :-0.0003304
Mean : 0.000000 Mean : 0.0000000 Mean : 0.0000000
3rd Qu.: 0.031472 3rd Qu.: 0.0263951 3rd Qu.: 0.0304928
Max. : 0.115554 Max. : 0.2985828 Max. : 0.2642922
vec11 vec19 vec5
Min. :-0.2190326 Min. :-0.1735843 Min. :-0.224152
1st Qu.:-0.0184901 1st Qu.:-0.0303940 1st Qu.:-0.026723
Median : 0.0007727 Median :-0.0002972 Median :-0.007429
Mean : 0.0000000 Mean : 0.0000000 Mean : 0.000000
3rd Qu.: 0.0151115 3rd Qu.: 0.0299778 3rd Qu.: 0.021318
Max. : 0.1915285 Max. : 0.1358331 Max. : 0.162426 sevQp <- fitted(meQp)
kreiseWithCovidMeldeWeeklyCleanedPredictorsAddedSvem <- st_sf(data.frame(kreiseWithCovidMeldeWeeklyCleanedPredictorsAdded,
sevQp))tm_shape(kreiseWithCovidMeldeWeeklyCleanedPredictorsAddedSvem) + tm_polygons(col = colnames(sevQp),
palette = "-RdBu", lwd = 0.5, n = 6, midpoint = 0, legend.show = FALSE) + tm_layout(main.title = "Selected spatial eigenvectors",
legend.outside = TRUE, attr.outside = TRUE, panel.show = TRUE, panel.labels = colnames(sevQp)) +
tm_scale_bar()
Rechecking spatial autocorrelation
(moranSevmQp <- lm.morantest(modelSevmQp, listw = unionedListw_d))
Global Moran I for regression residuals
data:
model: glm(formula = sumCasesTotal ~ Stimmenanteile.AfD +
Auslaenderanteil + Hochqualifizierte + fitted(meQp), family =
quasipoisson, data = kreiseWithCovidMeldeWeeklyCleanedPredictorsAdded,
offset = log(EWZ))
weights: unionedListw_d
Moran I statistic standard deviate = 7.3156, p-value = 1.281e-13
alternative hypothesis: greater
sample estimates:
Observed Moran I Expectation Variance
1.680228e-01 -2.875261e-06 5.275358e-04 We see that were is still some left over spatial autocorrelation not absorbed by the spatial eigenvectors. However, the amount of spatial autocorrelation has been reduced by a strong degree, from 0.51 to 0.17.
kreiseCentroids$residSevmQp <- residuals(modelSevmQp)
kreiseCentroids$residSevmQpAbs <- abs(kreiseCentroids$residSevmQp)
tm_shape(kreiseWithCovidMeldeWeeklyCleanedPredictorsAdded) + tm_polygons(col = "grey") +
tm_shape(kreiseCentroids) + tm_bubbles(size = "residSevmQpAbs", col = "residSevmQp",
palette = "-RdBu", alpha = 0.9, perceptual = TRUE, scale = 0.8, border.alpha = 0.3,
title.size = "Abs residual", title.col = "Residuals", n = 5) + tm_layout(main.title = "Pearson residuals, SEVM NB",
bg = "darkgrey", legend.outside = TRUE, attr.outside = TRUE) + tm_scale_bar()
## Spatial varying coefficients?
Which eigenvectors are correlated with ruralness? As ruralness became insignificant after incorporation of the spatial eigenvectors we might suspect that it is correlated with a subset of the eigenvectors. Let’s find out which eigenvectors are involved.
kreiseWithCovidMeldeWeeklyCleanedPredictorsAddedSvem %>%
st_drop_geometry() %>%
dplyr::select(Laendlichkeit, colnames(sevQp)) %>%
ggpairs()
Eigenvectors 2 and 1 are negatively correlated with ruralness (“Laendlichkeit”).
It might be that the spatial eigenvectors could be used to moderate the effect of ruralness - i.e. that the regression coefficient varies in space. Likewise we could investigate moderating relationships between other predictors and spatial eigenvectors.
colnames(sevQp) [1] "vec2" "vec6" "vec4" "vec10" "vec1" "vec18" "vec21" "vec11" "vec19"
[10] "vec5" modelSevmQpInt <- glm(sumCasesTotal ~ Stimmenanteile.AfD + Auslaenderanteil + Hochqualifizierte +
(vec1 + vec2) * Laendlichkeit + vec4 + vec6 + vec18 + vec21 + vec5 + vec11 +
vec19 + vec10, family = quasipoisson, offset = log(EWZ), data = kreiseWithCovidMeldeWeeklyCleanedPredictorsAddedSvem)
summary(modelSevmQpInt)
Call:
glm(formula = sumCasesTotal ~ Stimmenanteile.AfD + Auslaenderanteil +
Hochqualifizierte + (vec1 + vec2) * Laendlichkeit + vec4 +
vec6 + vec18 + vec21 + vec5 + vec11 + vec19 + vec10, family = quasipoisson,
data = kreiseWithCovidMeldeWeeklyCleanedPredictorsAddedSvem,
offset = log(EWZ))
Deviance Residuals:
Min 1Q Median 3Q Max
-56.457 -10.611 -1.923 7.039 64.584
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -3.7142500 0.0462340 -80.336 < 2e-16 ***
Stimmenanteile.AfD 0.0420443 0.0019723 21.317 < 2e-16 ***
Auslaenderanteil 0.0292675 0.0022308 13.120 < 2e-16 ***
Hochqualifizierte -0.0048261 0.0018138 -2.661 0.008125 **
vec1 -0.2840583 0.1984612 -1.431 0.153158
vec2 0.8584225 0.2093412 4.101 5.04e-05 ***
Laendlichkeit -0.0003932 0.0004360 -0.902 0.367675
vec4 -1.5274509 0.1610968 -9.482 < 2e-16 ***
vec6 -0.8291938 0.1780536 -4.657 4.43e-06 ***
vec18 -0.6106731 0.1490572 -4.097 5.11e-05 ***
vec21 -1.3629828 0.1711811 -7.962 1.94e-14 ***
vec5 0.5642140 0.1649767 3.420 0.000694 ***
vec11 -0.3582028 0.1479799 -2.421 0.015958 *
vec19 -0.9059447 0.1708088 -5.304 1.92e-07 ***
vec10 0.0975937 0.1615422 0.604 0.546110
vec1:Laendlichkeit -0.0167823 0.0076383 -2.197 0.028610 *
vec2:Laendlichkeit 0.0152535 0.0070868 2.152 0.031992 *
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
(Dispersion parameter for quasipoisson family taken to be 221.7727)
Null deviance: 344479 on 399 degrees of freedom
Residual deviance: 84440 on 383 degrees of freedom
AIC: NA
Number of Fisher Scoring iterations: 4drop1(modelSevmQpInt, test = "F")Single term deletions
Model:
sumCasesTotal ~ Stimmenanteile.AfD + Auslaenderanteil + Hochqualifizierte +
(vec1 + vec2) * Laendlichkeit + vec4 + vec6 + vec18 + vec21 +
vec5 + vec11 + vec19 + vec10
Df Deviance F value Pr(>F)
<none> 84440
Stimmenanteile.AfD 1 184596 454.2864 < 2.2e-16 ***
Auslaenderanteil 1 121710 169.0505 < 2.2e-16 ***
Hochqualifizierte 1 86012 7.1308 0.0079000 **
vec4 1 104417 90.6136 < 2.2e-16 ***
vec6 1 89229 21.7250 4.351e-06 ***
vec18 1 88181 16.9708 4.654e-05 ***
vec21 1 98821 65.2300 8.757e-15 ***
vec5 1 87016 11.6855 0.0006975 ***
vec11 1 85735 5.8741 0.0158275 *
vec19 1 90660 28.2149 1.845e-07 ***
vec10 1 84521 0.3672 0.5448907
vec1:Laendlichkeit 1 85511 4.8602 0.0280777 *
vec2:Laendlichkeit 1 85465 4.6505 0.0316661 *
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 11 - modelSevmQpInt$deviance/modelSevmQpInt$null.deviance[1] 0.7548771 - modelSevmQp$deviance/modelSevmQp$null.deviance[1] 0.7430323There is a significant interaction between ruralness and spatial eigenvector 2 and 1 while the main effect of ruralness is clearly not significant. However, the explained deviance does not increase much. How does the interaction look like in space?
mRuralVec2 <- kreiseWithCovidMeldeWeeklyCleanedPredictorsAddedSvem %>%
mutate(vec2Ruralness = vec2 * Laendlichkeit) %>%
tm_shape() + tm_polygons(col = c("vec2Ruralness"), legend.hist = TRUE, palette = "-RdBu",
style = "fisher", n = 6, midpoint = 0, legend.reverse = TRUE, title = "Ruralness moderated by eigenvector 2") +
tm_layout(legend.outside = TRUE, bg.color = "darkgrey", outer.bg.color = "lightgrey",
legend.outside.size = 0.5, attr.outside = TRUE, legend.outside.position = "left",
main.title = "Ruralness by sev 2") + tm_scale_bar()
tmap_arrange(mRuralVec2, mRural)
mRuralVec2 <- kreiseWithCovidMeldeWeeklyCleanedPredictorsAddedSvem %>%
mutate(vec2Ruralness = vec1 * Laendlichkeit) %>%
tm_shape() + tm_polygons(col = c("vec2Ruralness"), legend.hist = TRUE, palette = "-RdBu",
style = "fisher", n = 6, midpoint = 0, legend.reverse = TRUE, title = "Ruralness moderated by eigenvector 1") +
tm_layout(legend.outside = TRUE, bg.color = "darkgrey", outer.bg.color = "lightgrey",
legend.outside.size = 0.5, attr.outside = TRUE, legend.outside.position = "left",
main.title = "Ruralness by sev 1") + tm_scale_bar()
tmap_arrange(mRuralVec2, mRural)
lm.morantest(modelSevmQpInt, listw = unionedListw_d)
Global Moran I for regression residuals
data:
model: glm(formula = sumCasesTotal ~ Stimmenanteile.AfD +
Auslaenderanteil + Hochqualifizierte + (vec1 + vec2) * Laendlichkeit +
vec4 + vec6 + vec18 + vec21 + vec5 + vec11 + vec19 + vec10, family =
quasipoisson, data =
kreiseWithCovidMeldeWeeklyCleanedPredictorsAddedSvem, offset =
log(EWZ))
weights: unionedListw_d
Moran I statistic standard deviate = 6.9223, p-value = 2.222e-12
alternative hypothesis: greater
sample estimates:
Observed Moran I Expectation Variance
1.602309e-01 -2.989182e-06 5.358083e-04 kreiseCentroids$residSevmQpInt <- residuals(modelSevmQpInt)
kreiseCentroids$residSevmQpIntAbs <- abs(kreiseCentroids$residSevmQpInt)
tm_shape(kreiseWithCovidMeldeWeeklyCleanedPredictorsAdded) + tm_polygons(col = "grey") +
tm_shape(kreiseCentroids) + tm_bubbles(size = "residSevmQpIntAbs", col = "residSevmQpInt",
palette = "-RdBu", alpha = 0.9, perceptual = TRUE, scale = 0.8, border.alpha = 0.3,
title.size = "Abs residual", title.col = "Residuals", n = 5) + tm_layout(main.title = "Pearson residuals, SEVM QB\nruralness spatial varying coefficient",
bg = "darkgrey", legend.outside = TRUE, attr.outside = TRUE) + tm_scale_bar()
modelSevmQpInt2 <- update(modelSevmQpInt, ~. + Auslaenderanteil:vec2)
summary(modelSevmQpInt2)
Call:
glm(formula = sumCasesTotal ~ Stimmenanteile.AfD + Auslaenderanteil +
Hochqualifizierte + vec1 + vec2 + Laendlichkeit + vec4 +
vec6 + vec18 + vec21 + vec5 + vec11 + vec19 + vec10 + vec1:Laendlichkeit +
vec2:Laendlichkeit + Auslaenderanteil:vec2, family = quasipoisson,
data = kreiseWithCovidMeldeWeeklyCleanedPredictorsAddedSvem,
offset = log(EWZ))
Deviance Residuals:
Min 1Q Median 3Q Max
-53.822 -10.278 -1.923 7.356 58.824
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -3.823e+00 4.928e-02 -77.579 < 2e-16 ***
Stimmenanteile.AfD 4.603e-02 2.045e-03 22.512 < 2e-16 ***
Auslaenderanteil 3.303e-02 2.280e-03 14.490 < 2e-16 ***
Hochqualifizierte -3.396e-03 1.778e-03 -1.910 0.056931 .
vec1 -2.644e-01 1.919e-01 -1.377 0.169174
vec2 3.858e+00 5.901e-01 6.538 2.00e-10 ***
Laendlichkeit 4.514e-05 4.304e-04 0.105 0.916537
vec4 -1.513e+00 1.557e-01 -9.720 < 2e-16 ***
vec6 -7.843e-01 1.732e-01 -4.529 7.93e-06 ***
vec18 -6.628e-01 1.443e-01 -4.593 5.94e-06 ***
vec21 -1.208e+00 1.674e-01 -7.218 2.88e-12 ***
vec5 5.502e-01 1.597e-01 3.444 0.000636 ***
vec11 -5.739e-01 1.497e-01 -3.833 0.000148 ***
vec19 -8.453e-01 1.659e-01 -5.095 5.49e-07 ***
vec10 3.387e-01 1.628e-01 2.080 0.038173 *
vec1:Laendlichkeit -1.555e-02 7.412e-03 -2.097 0.036607 *
vec2:Laendlichkeit -1.486e-02 8.818e-03 -1.685 0.092810 .
Auslaenderanteil:vec2 -2.104e-01 3.875e-02 -5.428 1.02e-07 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
(Dispersion parameter for quasipoisson family taken to be 207.0252)
Null deviance: 344479 on 399 degrees of freedom
Residual deviance: 78251 on 382 degrees of freedom
AIC: NA
Number of Fisher Scoring iterations: 41 - modelSevmQpInt2$deviance/modelSevmQpInt2$null.deviance[1] 0.7728419If we include the interaction between share of foreigners and sev2 the interaction between ruralness and sev2 becomes insignificant.
modelSevmQpInt2 <- update(modelSevmQpInt2, ~. - vec2:Laendlichkeit)
summary(modelSevmQpInt2)
Call:
glm(formula = sumCasesTotal ~ Stimmenanteile.AfD + Auslaenderanteil +
Hochqualifizierte + vec1 + vec2 + Laendlichkeit + vec4 +
vec6 + vec18 + vec21 + vec5 + vec11 + vec19 + vec10 + vec1:Laendlichkeit +
Auslaenderanteil:vec2, family = quasipoisson, data = kreiseWithCovidMeldeWeeklyCleanedPredictorsAddedSvem,
offset = log(EWZ))
Deviance Residuals:
Min 1Q Median 3Q Max
-55.114 -10.601 -2.052 7.540 60.654
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -3.8167201 0.0491325 -77.682 < 2e-16 ***
Stimmenanteile.AfD 0.0461689 0.0020456 22.570 < 2e-16 ***
Auslaenderanteil 0.0323683 0.0022466 14.408 < 2e-16 ***
Hochqualifizierte -0.0035162 0.0017799 -1.975 0.048930 *
vec1 -0.2335194 0.1915538 -1.219 0.223564
vec2 3.1887116 0.4350967 7.329 1.39e-12 ***
Laendlichkeit 0.0001466 0.0004251 0.345 0.730472
vec4 -1.5321164 0.1556024 -9.846 < 2e-16 ***
vec6 -0.7651248 0.1733036 -4.415 1.32e-05 ***
vec18 -0.6187969 0.1422214 -4.351 1.74e-05 ***
vec21 -1.2094141 0.1679333 -7.202 3.17e-12 ***
vec5 0.5453670 0.1601063 3.406 0.000728 ***
vec11 -0.5547785 0.1493788 -3.714 0.000234 ***
vec19 -0.8707075 0.1655418 -5.260 2.40e-07 ***
vec10 0.2826184 0.1597694 1.769 0.077704 .
vec1:Laendlichkeit -0.0116193 0.0070251 -1.654 0.098950 .
Auslaenderanteil:vec2 -0.1693540 0.0301058 -5.625 3.58e-08 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
(Dispersion parameter for quasipoisson family taken to be 207.8777)
Null deviance: 344479 on 399 degrees of freedom
Residual deviance: 78840 on 383 degrees of freedom
AIC: NA
Number of Fisher Scoring iterations: 41 - modelSevmQpInt2$deviance/modelSevmQpInt2$null.deviance[1] 0.7711321modelSevmQpInt2 <- update(modelSevmQpInt2, ~. - vec1:Laendlichkeit)
summary(modelSevmQpInt2)
Call:
glm(formula = sumCasesTotal ~ Stimmenanteile.AfD + Auslaenderanteil +
Hochqualifizierte + vec1 + vec2 + Laendlichkeit + vec4 +
vec6 + vec18 + vec21 + vec5 + vec11 + vec19 + vec10 + Auslaenderanteil:vec2,
family = quasipoisson, data = kreiseWithCovidMeldeWeeklyCleanedPredictorsAddedSvem,
offset = log(EWZ))
Deviance Residuals:
Min 1Q Median 3Q Max
-55.008 -10.009 -1.998 7.461 58.275
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -3.8276546 0.0487983 -78.438 < 2e-16 ***
Stimmenanteile.AfD 0.0469824 0.0019884 23.629 < 2e-16 ***
Auslaenderanteil 0.0327520 0.0022406 14.617 < 2e-16 ***
Hochqualifizierte -0.0036675 0.0017805 -2.060 0.040093 *
vec1 -0.4029965 0.1614259 -2.496 0.012962 *
vec2 3.4764206 0.3998051 8.695 < 2e-16 ***
Laendlichkeit 0.0003002 0.0004144 0.724 0.469226
vec4 -1.6344266 0.1428255 -11.444 < 2e-16 ***
vec6 -0.7134649 0.1713180 -4.165 3.86e-05 ***
vec18 -0.6291195 0.1424453 -4.417 1.31e-05 ***
vec21 -1.2150457 0.1681284 -7.227 2.69e-12 ***
vec5 0.5278639 0.1600945 3.297 0.001068 **
vec11 -0.5757764 0.1492522 -3.858 0.000134 ***
vec19 -0.9210587 0.1630916 -5.647 3.17e-08 ***
vec10 0.3039040 0.1595634 1.905 0.057579 .
Auslaenderanteil:vec2 -0.1847810 0.0287029 -6.438 3.62e-10 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
(Dispersion parameter for quasipoisson family taken to be 208.6265)
Null deviance: 344479 on 399 degrees of freedom
Residual deviance: 79408 on 384 degrees of freedom
AIC: NA
Number of Fisher Scoring iterations: 41 - modelSevmQpInt2$deviance/modelSevmQpInt2$null.deviance[1] 0.7694841modelSevmQpInt2 <- update(modelSevmQpInt2, ~. - Laendlichkeit)
summary(modelSevmQpInt2)
Call:
glm(formula = sumCasesTotal ~ Stimmenanteile.AfD + Auslaenderanteil +
Hochqualifizierte + vec1 + vec2 + vec4 + vec6 + vec18 + vec21 +
vec5 + vec11 + vec19 + vec10 + Auslaenderanteil:vec2, family = quasipoisson,
data = kreiseWithCovidMeldeWeeklyCleanedPredictorsAddedSvem,
offset = log(EWZ))
Deviance Residuals:
Min 1Q Median 3Q Max
-55.793 -9.913 -2.073 7.366 59.329
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -3.805067 0.037489 -101.498 < 2e-16 ***
Stimmenanteile.AfD 0.046717 0.001954 23.908 < 2e-16 ***
Auslaenderanteil 0.032204 0.002108 15.277 < 2e-16 ***
Hochqualifizierte -0.004154 0.001648 -2.521 0.012113 *
vec1 -0.456773 0.143218 -3.189 0.001543 **
vec2 3.384027 0.378769 8.934 < 2e-16 ***
vec4 -1.639446 0.142565 -11.500 < 2e-16 ***
vec6 -0.712172 0.171235 -4.159 3.94e-05 ***
vec18 -0.625980 0.142251 -4.401 1.40e-05 ***
vec21 -1.222954 0.167732 -7.291 1.76e-12 ***
vec5 0.525974 0.159953 3.288 0.001101 **
vec11 -0.559127 0.147342 -3.795 0.000172 ***
vec19 -0.921057 0.163058 -5.649 3.15e-08 ***
vec10 0.294258 0.158882 1.852 0.064784 .
Auslaenderanteil:vec2 -0.178247 0.027241 -6.543 1.92e-10 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
(Dispersion parameter for quasipoisson family taken to be 208.4964)
Null deviance: 344479 on 399 degrees of freedom
Residual deviance: 79517 on 385 degrees of freedom
AIC: NA
Number of Fisher Scoring iterations: 41 - modelSevmQpInt2$deviance/modelSevmQpInt2$null.deviance[1] 0.7691665m1 <- kreiseWithCovidMeldeWeeklyCleanedPredictorsAddedSvem %>%
mutate(vec2foreigners = vec2 * Auslaenderanteil) %>%
tm_shape() + tm_polygons(col = c("vec2foreigners"), legend.hist = TRUE, palette = "-RdBu",
style = "fisher", n = 6, midpoint = 0, legend.reverse = TRUE, title = "Foreigners moderated by eigenvector 2") +
tm_layout(legend.outside = TRUE, bg.color = "darkgrey", outer.bg.color = "lightgrey",
legend.outside.size = 0.5, attr.outside = TRUE, legend.outside.position = "left",
main.title = "Foreigners by sev 2") + tm_scale_bar()
tmap_arrange(m1, mForeigner)
The effect of the share of foreigners differs in space. A potential explanation would be that the effect of the foreign population is moderated by education (at the individual level) as well as by language skills and economic status (of course all factors are interlinked). Also the measures to reach non german speaking inhabitants differs presumably between federal states - leading to potential differences in how those groups obey to rules such as mask wearing, social distancing or with respect to vaccination rates.
Of course there is much more potential to explore relationships between predictors and spatial eigenvectors and we have not even touched on interactions between “normal” predictors or higher order effects. We could (and should) try other predictor eigenvector combinations as well as interactions between the predictors.